Spectroscopic Data Reduction

Spectroscopic data reduction is the process of converting raw spectroscopic images from a spectrograph into a calibrated spectrum. This process involves several steps, including bias subtraction, flat fielding, wavelength calibration, and extraction of the spectrum. In this notebook, we will go through the process of reducing a spectroscopic image using Python.

The CCD pre-processing steps are as usual: bias/dark subtraction, flat fielding, and cosmic ray removal. The spectroscopic data reduction steps are: extraction of the spectrum, wavelength calibration and flux calibration. Here we will focus mostly on the spectroscopic data reduction steps. After the data reduction, the spectrum can be used for scientific analysis, such as measuring redshifts, line identifications, and measuring line fluxes. In this notebook, we will reduce a 2D spectrum of an afterglow of a gamma-ray burst (GRB). Since this tutorial is prepared for the Project No. 5 of our TAO class, I leave the scientific analysis of the spectrum to readers.

The data used in this notebook is from the Subaru Telescope’s Faint Object Camera and Spectrograph (FOCAS). The data is 2D images of a spectrum of an afterglow of a GRB, GRB 130606A, which is a significantly redshifted source at \(z\sim5.9\). The data was taken on June 7, 2013, using the VPH900 grism, which covers the wavelength range of 7500-10450 Angstroms (Subaru FOCAS Hompage).

Let’s start by organizing the data and displaying the 2D spectral images.

import warnings
from datetime import datetime
from pathlib import Path

import numpy as np
import preproc  # custom module for data reduction
from astropy.io import ascii, fits
from astropy.modeling.fitting import LevMarLSQFitter
from astropy.modeling.models import Chebyshev2D, Gaussian1D
from astropy.nddata import CCDData
from astropy.stats import gaussian_fwhm_to_sigma, sigma_clip
from astropy.table import Column, Table
from astropy.visualization import ZScaleInterval
from ccdproc import combine, cosmicray_lacosmic, gain_correct
from IPython.display import Image
from matplotlib import gridspec
from matplotlib import pyplot as plt
from matplotlib import rc, rcParams
from mpl_toolkits.axes_grid1 import make_axes_locatable
from numpy.polynomial.chebyshev import chebfit, chebval
from scipy.interpolate import UnivariateSpline
from scipy.ndimage import median_filter
from scipy.optimize import curve_fit
from scipy.signal import peak_widths
from skimage.feature import peak_local_max
from astropy import units as u

# plt.rcParams['axes.linewidth'] = 2
plt.rcParams["font.size"] = 12
plt.rcParams["axes.labelsize"] = 14

WD = Path("./data/proj5")
RAWDIR = WD / "raw" / "P03hbahk0304145443GZ_data"
OUTDIR = WD / "out"

if not OUTDIR.exists():
    OUTDIR.mkdir(parents=True)

warnings.filterwarnings("ignore", category=UserWarning, append=True)

def make_summary_table_focas(rawdir, suffix=".fits.gz"):
    """Make a summary table of the raw data in the directory. This function is
    specifically designed for the data taken with the FOCAS instrument.

    Args:
        rawdir (pathlib.Path): The directory containing the raw data.
        suffix (str, optional): The suffix of the raw data files. Defaults to
            '.fits.gz'.

    Returns:
        stab (astropy.table.Table): The summary table of the raw data in the directory.

    """
    # making a summary table
    summary = []
    for f in rawdir.glob("*" + suffix):
        hdr = fits.getheader(f)

        # getting the filter information from the header
        filt1 = hdr["FILTER01"]
        filt2 = hdr["FILTER02"]
        filt3 = hdr["FILTER03"]

        X = hdr["AIRMASS"]  # airmass
        disp = hdr["DISPERSR"]  # disperser

        # appending the data to the summary list
        summary.append(
            [
                f.name,
                hdr["DATE-OBS"],
                hdr["OBJECT"],
                hdr["RA"],
                hdr["DEC"],
                hdr["DATA-TYP"],
                hdr["EXPTIME"],
                X,
                hdr["ALTITUDE"],
                disp,
                hdr["SLT-WID"],
                hdr["SLT-PA"],
                hdr["SLIT"],
                filt1,
                filt2,
                filt3,
                hdr["UT-STR"],
                hdr["UT-END"],
            ]
        )

    # creating the summary table
    stab = Table(
        rows=summary,
        names=[
            "filename",
            "date-obs",
            "object",
            "ra",
            "dec",
            "type",
            "exptime",
            "airmass",
            "altitude",
            "disperser",
            "slit-width",
            "slit-pa",
            "slit",
            "filter1",
            "filter2",
            "filter3",
            "ut-str",
            "ut-end",
        ],
        dtype=[
            "U50",
            "U50",
            "U50",
            "U50",
            "U50",
            "U50",
            "f8",
            "f8",
            "f8",
            "U50",
            "f8",
            "f8",
            "U50",
            "U50",
            "U50",
            "U50",
            "U50",
            "U50",
        ],
    )
    return stab

summary_table = make_summary_table_focas(RAWDIR)
summary_table

Table length=60

filename	date-obs	object	ra	dec	type	exptime	airmass	altitude	disperser	slit-width	slit-pa	slit	filter1	filter2	filter3	ut-str	ut-end
str50	str50	str50	str50	str50	str50	float64	float64	float64	str50	float64	float64	str50	str50	str50	str50	str50	str50
FCSA00141794.fits.gz	2013-06-07	DOMEFLAT	11:17:48.038	+19:54:08.26	DOMEFLAT	3.0	1.0	89.9361	SCFCGRHD90	0.5	-90.3	SCFCSLLC08	NONE	SCFCFLSO58	NONE	04:37:11.906	04:37:15.139
FCSA00141920.fits.gz	2013-06-07	GRB 130606A	16:37:41.157	+29:47:45.71	OBJECT	1200.0	1.136	61.6129	SCFCGRHD90	0.5	-90.3	SCFCSLLC08	NONE	SCFCFLSO58	NONE	07:59:00.676	08:19:01.117
FCSA00141930.fits.gz	2013-06-07	GRB 130606A	16:37:41.230	+29:47:45.13	OBJECT	1200.0	1.038	74.4191	SCFCGRHD90	0.5	-90.3	SCFCSLLC08	NONE	SCFCFLSO58	NONE	09:03:28.980	09:23:29.405
FCSA00141848.fits.gz	2013-06-07	BIAS	11:39:50.841	+19:54:11.18	BIAS	0.0	1.0	89.9361	SCFCGRMR01	0.5	-90.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	04:59:10.571	04:59:10.574
FCSA00141916.fits.gz	2013-06-07	GRB 130606A	16:37:41.020	+29:47:44.75	OBJECT	1200.0	1.254	52.8249	SCFCGRHD90	0.5	-90.3	SCFCSLLC08	NONE	SCFCFLSO58	NONE	07:18:00.165	07:38:00.608
FCSA00141734.fits.gz	2013-06-06	BIAS	21:41:03.767	+17:46:28.41	BIAS	0.0	1.001	86.8997	SCFCGRHD90	0.5	-0.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	14:53:03.807	14:53:03.810
FCSA00141566.fits.gz	2013-06-06	DOMEFLAT	11:23:03.374	+19:54:09.25	DOMEFLAT	1.1	1.0	89.9625	SCFCGRHD90	0.5	-90.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	04:46:28.971	04:46:30.721
FCSA00141576.fits.gz	2013-06-06	DOMEFLAT	11:26:18.779	+19:54:09.76	DOMEFLAT	3.0	1.0	89.9625	SCFCGRHD90	0.5	-90.3	SCFCSLLC08	NONE	SCFCFLSO58	NONE	04:49:43.639	04:49:46.993
FCSA00141822.fits.gz	2013-06-07	DOMEFLAT	11:26:05.681	+19:54:09.65	DOMEFLAT	1.1	1.0	89.9361	SCFCGRHD90	0.5	-90.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	04:45:27.990	04:45:29.664
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
FCSA00141730.fits.gz	2013-06-06	BIAS	23:03:44.623	-01:20:01.29	BIAS	0.0	1.169	58.8068	SCFCGRHD90	0.5	-0.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	14:52:13.406	14:52:13.410
FCSA00141818.fits.gz	2013-06-07	DOMEFLAT	11:25:08.002	+19:54:09.51	DOMEFLAT	1.1	1.0	89.9361	SCFCGRHD90	0.5	-90.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	04:44:30.393	04:44:31.697
FCSA00142010.fits.gz	2013-06-07	Feige110	23:20:01.964	-05:10:51.68	OBJECT	120.0	1.208	55.8712	SCFCGRHD90	0.5	-405.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	15:03:10.694	15:05:11.030
FCSA00141798.fits.gz	2013-06-07	DOMEFLAT	11:18:59.560	+19:54:08.48	DOMEFLAT	3.0	1.0	89.9361	SCFCGRHD90	0.5	-90.3	SCFCSLLC08	NONE	SCFCFLSO58	NONE	04:38:23.184	04:38:26.574
FCSA00141854.fits.gz	2013-06-07	BIAS	11:41:06.075	+19:54:11.27	BIAS	0.0	1.0	89.9361	SCFCGRMR01	0.5	-90.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	05:00:25.637	05:00:25.674
FCSA00141586.fits.gz	2013-06-06	BIAS	11:32:20.002	+19:54:10.54	BIAS	0.0	1.0	89.9625	SCFCGRHD90	0.5	-90.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	04:55:43.880	04:55:43.883
FCSA00141562.fits.gz	2013-06-06	DOMEFLAT	11:22:06.397	+19:54:09.10	DOMEFLAT	1.1	1.0	89.9625	SCFCGRHD90	0.5	-90.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	04:45:32.080	04:45:33.717
FCSA00141572.fits.gz	2013-06-06	DOMEFLAT	11:25:16.188	+19:54:09.60	DOMEFLAT	3.0	1.0	89.9625	SCFCGRHD90	0.5	-90.3	SCFCSLLC08	NONE	SCFCFLSO58	NONE	04:48:41.382	04:48:44.710
FCSA00141728.fits.gz	2013-06-06	BIAS	23:19:58.381	-05:11:14.89	BIAS	0.0	1.248	53.1982	SCFCGRHD90	0.5	-0.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	14:51:48.193	14:51:48.197
FCSA00141800.fits.gz	2013-06-07	DOMEFLAT	11:19:30.656	+19:54:08.58	DOMEFLAT	3.0	1.0	89.9361	SCFCGRHD90	0.5	-90.3	SCFCSLLC08	NONE	SCFCFLSO58	NONE	04:38:54.202	04:38:57.612

biastab = summary_table[summary_table["type"] == "BIAS"]
flattab = summary_table[summary_table["type"] == "DOMEFLAT"]
arctab = summary_table[
    (summary_table["type"] == "COMPARISON") & (summary_table["slit"] == "SCFCSLLC08")
]
scitab = summary_table[summary_table["object"] == "GRB 130606A"]
stdtab_f110 = summary_table[
    (summary_table["object"] == "Feige110") & (summary_table["type"] == "OBJECT")
]
stdtab_f34 = summary_table[
    (summary_table["object"] == "Feige34") & (summary_table["type"] == "OBJECT")
]

# for displaying in zscale
interval = ZScaleInterval()

# plotting the images
fig, axes = plt.subplots(2, 2, figsize=(15, 4))
titles = ["Bias", "Flat", "Arc", "Sci"]
bias = fits.getdata(RAWDIR / biastab["filename"][0])
flat = fits.getdata(RAWDIR / flattab["filename"][0])
arc = fits.getdata(RAWDIR / arctab["filename"][0])
sci = fits.getdata(RAWDIR / scitab["filename"][0])
for i, img in enumerate([bias, flat, arc, sci]):
    vmin, vmax = interval.get_limits(img)
    ax = axes[i // 2][i % 2]
    ax.imshow(img.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
    ax.text(
        0.05,
        0.95,
        titles[i],
        transform=ax.transAxes,
        ha="left",
        va="top",
        bbox=dict(facecolor="white", alpha=0.8),
    )
    ax.tick_params(axis="both", length=0.0, labelleft=False, labelbottom=False)
plt.tight_layout()

del bias, flat, arc, sci

Preprocessing of the Raw Images

Preprocessing of the raw images is the first step in the data reduction process. The raw images are usually not ready for scientific analysis and need to be preprocessed before they can be used for scientific analysis. The preprocessing steps include the following:

1. Bias subtraction

2. Dark subtraction

3. Flat fielding

4. Cosmic ray removal

There are lots of software available for preprocessing the raw images. In this tutorial, we will use the ccdproc package for preprocessing the raw images. The ccdproc package is a Python package for reducing astronomical data taken with CCD cameras. It provides many of the necessary tools for processing of CCD images. For more information on the ccdproc package, please visit the documentation, and for the complete guide on the reduction process for astronomical images using Python, CCD Reduction and Photometry Guide.

In this notebook, we will preprocess the raw images using preproc.py script. This script is available in the tutorial directory. Dark subtraction will not be performed in this notebook since the amount of dark current in the images is negligible.

Making the Master Bias

Before making the master bias, we need to check the bias images for any defects. We skip this step here, but it is important to check since the bias images are used to calibrate all the other images. If there are defects in the bias images, they will be propagated to all the other images.

The master bias is made by combining all the bias images. The master bias will be subtracted from all the images to remove the bias level. For combining the bias images, the sigma-clipping algorithm is used to reject outliers and then the average of the remaining images is calculated.

biaslist = [RAWDIR / f for f in biastab["filename"]]
_ = preproc.combine_bias(biaslist, outdir=OUTDIR, outname="MBias.fits")
del _

Show code cell output Hide code cell output

INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]

WARNING: FITSFixedWarning: PC001001= 1.00000000 / Pixel Coordinate translation matrix 
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC001002= 0.00000000 / Pixel Coordinate translation matrix 
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC002001= 0.00000000 / Pixel Coordinate translation matrix 
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC002002= 1.00000000 / Pixel Coordinate translation matrix 
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / The equatorial coordinate system 
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T04:59:10.645' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'unitfix' made the change 'Changed units:
  'degree' -> 'deg'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'celfix' made the change 'Unmatched celestial axes'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T14:53:03.880' from MJD-END'. [astropy.wcs.wcs]

INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]

WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:56:58.652' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T04:59:35.780' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T15:06:35.545' from MJD-END'. [astropy.wcs.wcs]

INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]

WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T05:00:00.751' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:56:34.054' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T15:07:00.686' from MJD-END'. [astropy.wcs.wcs]

INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]

WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T15:07:53.456' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:56:09.042' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T14:53:29.035' from MJD-END'. [astropy.wcs.wcs]

INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]

WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T15:06:10.558' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:55:19.314' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T04:58:45.726' from MJD-END'. [astropy.wcs.wcs]

INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]

WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T14:52:38.925' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T15:07:26.756' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T14:52:13.481' from MJD-END'. [astropy.wcs.wcs]

INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]

WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T05:00:25.744' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:55:43.954' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T14:51:48.267' from MJD-END'. [astropy.wcs.wcs]

INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]

mbias = OUTDIR / "MBias.fits"
mbias_img = fits.getdata(mbias)
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
vmin, vmax = interval.get_limits(mbias_img)
im = ax.imshow(mbias_img.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Master Bias")
divider = make_axes_locatable(ax)
cax = divider.append_axes("bottom", size="5%", pad=0.5)
plt.colorbar(im, cax=cax, orientation="horizontal")
fig.show()

Making Master Flat

Not like the photometric observation, spectroscopic flat image is typically taken with something like halogen lamp or tungsten lamp, which is inevitably not uniform especially for the dispersion direction. Therefore, we need to make a master flat image flat, which involves modeling the flat spectrum and dividing the flat image by the model. To obtain a good master flat, we need to:

1. Examine the flat spectra: do they show any variablity throughout the observation? can we combine them without any further consideration? do they have any features like wiggles or bumps? do they have any discontinuity?

2. Model the flat spectra: we can use a polynomial (or smoothing convolution) to model the flat spectrum. The degree of the polynomial (or size of the smoothing kernel) should be determined by the complexity of the flat spectrum.

3. Normalize the flat spectra: divide the flat image by the model to obtain a normalized flat image.

Combining can be done either before or after the modeling. If the flat spectra are stable, we can combine them first and then model the combined flat spectrum. If the flat spectra are not stable, we can model each flat spectrum and then combine the normalized flat images. If the flat spectra show significant variability, we can also consider not combining them at all.

Because of this complexity, making a master flat can be a time-consuming part of the data reduction process. Moreover, the quality of the master flat directly affects the quality of the reduced science images, which may introduce systematic errors if not done properly. Therefore, the spectroscopic flat fielding is often omitted in the data reduction process, when the flat fielding conversely introduces more errors than it corrects.

Checking for Flat Image Shift

The flat image may shift slightly depending on the telescope’s pointing direction. Since different images might require different flats, it can be risky to median combine all the flat images taken overnight. We should perform a quick and rough check to see if there is a shift in the flat images.

To avoid bias might arise from this shift, it is recommended to take flat images at the same pointing direction as the science images. However, if the shift is small, we can still use the median combined flat image for calibration.

Note

With the same reason, we should be careful when we reduce the arc frames. The location of the lines in the arc spectrum may shift slightly depending on the telescope’s pointing direction. Therefore, it is strongly recommended to take arc frames at the same pointing direction as the science frames.

We should check for the shift in the arc frames, but we will skip this process for simplicity.

However, unlike the flat fielding, the wavelength calibration is of critical importance in the spectroscopic data reduction process, since few percent of errors in the wavelength calibration can introduce systematic errors in the scientific analysis, such as the measurement of the redshift of galaxies. Therefore, we should always check for the shift in the arc frames before wavelength calibration.

In this notebook, we will skip this process for simplicity, but you should always check for flat image shift in your own data. We will assume that the flat images are already aligned, and we will proceed to combine them.

flattab_exp3 = flattab[flattab["exptime"] == 3.0]
flattab_exp1 = flattab[flattab["exptime"] == 1.1]
flatlist_exp3 = [RAWDIR / f for f in flattab_exp3["filename"]]
flatlist_exp1 = [RAWDIR / f for f in flattab_exp1["filename"]]

mbias = OUTDIR / "MBias.fits"
_ = preproc.make_master_flat(
    flatlist_exp3,
    outdir=OUTDIR,
    outname="MFlat03.0.fits",
    mbias=mbias,
    filter_key="FILTER02",
    verbose=False,
)
_ = preproc.make_master_flat(
    flatlist_exp1,
    outdir=OUTDIR,
    outname="MFlat01.1.fits",
    mbias=mbias,
    filter_key="FILTER02",
    verbose=False,
)
del _

Comparing with the Standard Star Profile

Typical flat spectra shows many wiggles and bumps. Yet, we don’t know if these features are originated from the grating (as a representation of common instruments that can affect the both flat and science frames), which in turn will affect our science frame equally. If they are not due to the grating, such as filters in front of the flat-field lamp, the extreme difference in color temperature between the lamps and the object of interest, or wavelength-dependent variations in the reflectivity of the paint used for the flat-field screen (Massey & Hansen, 2010), they are not going to present on the science frame.

Let’s see if a supposedly smooth source (in this case, a standard star) has same “wiggles” and “bumps” on its spectrum. If this is the case, we should normalize the flat spectrum with a constant value (or a “low-order” functional model). If not, we should adopt “high-order” model or constant to normalize our flat field. For more information, please refer to Section 3.1., Massey & Hansen (2010).

fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
flat03 = fits.getdata(OUTDIR / "MFlat03.0.fits")
vmin, vmax = interval.get_limits(flat03)
im = ax.imshow(flat03.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Master Flat (3.0s)")

fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
flat01 = fits.getdata(OUTDIR / "MFlat01.1.fits")
vmin, vmax = interval.get_limits(flat01)
im = ax.imshow(flat03.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Master Flat (1.1s)")

fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
f34 = fits.getdata(RAWDIR / stdtab_f34["filename"][0])
vmin, vmax = interval.get_limits(f34)
im = ax.imshow(f34.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Feige34 (Standard Star)")

Text(0.5, 1.0, 'Feige34 (Standard Star)')

fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(211)
ax.plot(flat03[:, 700])
ax.text(
    0.05,
    0.95,
    "Master Flat, 3.0s, Row 700",
    transform=ax.transAxes,
    ha="left",
    va="top",
)
ax.grid()
ax = fig.add_subplot(212)
ax.plot(f34[:, 737])
ax.text(
    0.05, 0.95, "Feige34, Row 700, Raw", transform=ax.transAxes, ha="left", va="top"
)
ax.grid()

Although we can see some tiny wiggles in the flat spectrum, this is not significant enough to adopt higher-order polynomials to model the flat spectrum, as opposed to the case if we have significant fringe patterns. If you closely examine the flat spectrum, you will see that the wiggles are not random, but rather located at the atmospheric absorption lines. This is because the flat image is taken with a lamp, and the lamp light is absorbed by the atmosphere at these absorption lines. Since the absorption lines are not random, we can safely ignore them when modeling the flat spectrum.

Here we model the flat spectrum with median filtering, which is a simple and effective way to remove noise while preserving the shape of the spectrum. We will use a median filter in the dispersion direction to model the flat spectrum, with a window size of 100 pixels. This is a good starting point, but you may need to adjust the window size depending on the complexity of the flat spectrum.

Note

Although I compared the horizontal profile of the standard star and the flat image for convenience, it might as well natural to compare the final extracted spectra of the standard star after applying both the high-order and low-order model of the flat field.

For simplicity, we will only use the combined flat image taken with the exposure time of 3 seconds. In practice, you may use all the flat images taken with different exposure times to make a master flat, or you can use the flat images taken with long and short exposure times to make a bad pixel map.

# median filtering the flat field image along the dispersion axis
flat03_filt = median_filter(flat03, size=50, axes=0)
# normalizing the flat field image
normflat = flat03 / flat03_filt

fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(211)
ax.plot(flat03[:, 700], c="k", label="Raw", lw=1)
ax.plot(flat03_filt[:, 700], c="r", label="Filtered", lw=0.8)
ax.set_xlim(500, 4100)
ax.set_ylim(0, 2.2)
ax.legend()
ax.grid()
ax = fig.add_subplot(212)
ax.plot(normflat[:, 700], c="k", label="Normalized", lw=1)
ax.set_xlim(500, 4100)
ax.set_ylim(0.95, 1.05)
ax.legend()
ax.grid()

fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(311)
vmin, vmax = interval.get_limits(flat03)
im = ax.imshow(flat03.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Master Flat")
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="2%", pad=0.05)
plt.colorbar(im, cax=cax)
ax = fig.add_subplot(312)
im = ax.imshow(flat03_filt.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Filtered Master Flat")
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="2%", pad=0.05)
plt.colorbar(im, cax=cax)
ax = fig.add_subplot(313)
vmin, vmax = interval.get_limits(normflat)
im = ax.imshow(normflat.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Normalized Master Flat")
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="2%", pad=0.05)
plt.colorbar(im, cax=cax)
fig.show()

# writing the normalized flat field image to a file
hdr = fits.getheader(OUTDIR / "MFlat03.0.fits")
hdu = fits.PrimaryHDU(normflat, header=hdr)
hdu.writeto(OUTDIR / "nMFlat03.0.fits", overwrite=True)

Preprocessing the Raw Images

Now we will preprocess the raw images. We will subtract the master bias and divide the master flat from the raw science frames. These preprocessed science frames will not be combined, since the spatial location of the object is different in each frame. We will combine the light from the object after extracting the 1D spectrum.

mbias = OUTDIR / "MBias.fits"
mflat = OUTDIR / "nMFlat03.0.fits"
# preprocessing the arc image
arc_list = [RAWDIR / f for f in arctab["filename"]]
print("Preprocessing arc images...")
arc = preproc.preproc(arc_list, mbias=mbias, mflat=mflat, outdir=OUTDIR)

# combine the preprocessed arc image
print("Combining the preprocessed arc images...")
parcs = [CCDData.read(OUTDIR / ("p" + f.stem)) for f in arc_list]
carc = combine(
    parcs,
    method="average",
    sigma_clip=True,
    sigma_clip_low_thresh=3,
    sigma_clip_high_thresh=3,
)

carc.write(OUTDIR / "cArc.fits", overwrite=True)
print("The combined arc image is saved as cArc.fits")

# preprocessing the science image
print("Preprocessing science images...")
sci_list = [RAWDIR / f for f in scitab["filename"]]
_ = preproc.preproc(sci_list, mbias=mbias, mflat=mflat, outdir=OUTDIR, insert_ivar=True)

# preprocessing the standard star images
print("Preprocessing standard star images...")
print("Standard star: Feige110")
std_list_f110 = [RAWDIR / f for f in stdtab_f110["filename"]]
_ = preproc.preproc(
    std_list_f110, mbias=mbias, mflat=mflat, outdir=OUTDIR, insert_ivar=True
)
print("Standard star: Feige34")
std_list_f34 = [RAWDIR / f for f in stdtab_f34["filename"]]
_ = preproc.preproc(
    std_list_f34, mbias=mbias, mflat=mflat, outdir=OUTDIR, insert_ivar=True
)

Show code cell output Hide code cell output

Preprocessing arc images...
Done: pFCSA00141840.fits

/Users/hbahk/class/TAOtutorials/tutorials/preproc.py:295: RuntimeWarning: invalid value encountered in divide
  sci_data /= mflat_data    # Flat fielding

Done: pFCSA00141882.fits
Done: pFCSA00141580.fits
Combining the preprocessed arc images...
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]

WARNING: FITSFixedWarning: PC001001= 1.00000000 / Pixel Coordinate translation matrix 
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC001002= 0.00000000 / Pixel Coordinate translation matrix 
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC002001= 0.00000000 / Pixel Coordinate translation matrix 
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC002002= 1.00000000 / Pixel Coordinate translation matrix 
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / The equatorial coordinate system 
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T04:55:27.734' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'unitfix' made the change 'Changed units:
  'degree' -> 'deg'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'celfix' made the change 'Unmatched celestial axes'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T05:51:12.249' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:52:48.140' from MJD-END'. [astropy.wcs.wcs]

INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
 [astropy.nddata.ccddata]
The combined arc image is saved as cArc.fits
Preprocessing science images...
Done: pFCSA00141920.fits

/Users/hbahk/class/TAOtutorials/tutorials/preproc.py:295: RuntimeWarning: invalid value encountered in divide
  sci_data /= mflat_data    # Flat fielding

Done: pFCSA00141930.fits
Done: pFCSA00141916.fits
Done: pFCSA00141928.fits
Done: pFCSA00141938.fits
Done: pFCSA00141940.fits
Done: pFCSA00141918.fits
Done: pFCSA00141936.fits
Done: pFCSA00141926.fits
Preprocessing standard star images...
Standard star: Feige110
Done: pFCSA00142004.fits
Done: pFCSA00142006.fits
Done: pFCSA00142010.fits
Standard star: Feige34
Done: pFCSA00141880.fits
Done: pFCSA00141894.fits
Done: pFCSA00141878.fits

fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
vmin, vmax = interval.get_limits(carc.data)
im = ax.imshow(carc.data.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Combined Image")

fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
sci = fits.getdata(OUTDIR / f"p{sci_list[0].stem}")
vmin, vmax = interval.get_limits(sci)
im = ax.imshow(sci.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Preprocessed Science Image")

fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
std34 = fits.getdata(OUTDIR / f"p{std_list_f34[0].stem}")
vmin, vmax = interval.get_limits(std34)
im = ax.imshow(std34.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Preprocessed Feige34 Image")
fig.show()

Extracting 1D Spectra

The 1D spectra are extracted from the 2D science frames by summing the flux along the spatial direction. The extraction process consists of the following steps:

1. Source Identification: Identify the source in the 2D science frame.

2. Aperture Tracing: Trace the spatial position of the source along the dispersion direction.

3. Aperture Summation: Sum the flux along the trace to extract the 1D spectrum.

Let’s extract the 1D spectra of the standard star first, because it is brighter and easier to extract.

1. Source Identification

Let’s trim the images to the region of interest. Then, after the cosmic ray removal, we will identify the source in the 2D science frame. We will use the peak_local_max function from skimage.feature module to identify the source in a slice of the 2D science frame. The peak_local_max function finds the local maxima, which can be used to identify the source in the 2D science frame.

# Trim the images
TRIM = [0, 4220, 670, 810]  # y_lower, y_upper, x_lower, x_upper
fname = OUTDIR / f"p{std_list_f34[0].stem}"
stdhdu = fits.open(fname)
stdhdr = stdhdu[0].header
stdtrim = stdhdu[0].data[TRIM[0] : TRIM[1], TRIM[2] : TRIM[3]]
exptime = stdhdr["EXPTIME"]
rdnoise = stdhdr["RDNOISE"]

fig = plt.figure(figsize=(12, 2))
ax = fig.add_subplot(111)
vmin, vmax = interval.get_limits(stdtrim)
im = ax.imshow(stdtrim.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Preprocessed Feige34 Image (Trimmed)")
fig.show()

The TRIM parameter is used to trim the images to the region of interest. The elements of the TRIM parameter is determined by examining the 2D science frame, to make the source centered in the image.

Cosmic ray removal often takes a long time. Hence, we will apply cosmic ray removal on the trimmed images to save unnecessary computation. To get reliable removal of cosmic rays, we need to mask bad pixels, but here we will skip this step for simplicity. For more information on cosmic ray removal, please refer to the CCD Data Reduction Guide.

# Cosmic ray removal
stdccd = CCDData(data=stdtrim, header=stdhdr, unit="adu")
gcorr = gain_correct(stdccd, gain=stdhdr["GAIN"] * u.electron / u.adu)
crrej = cosmicray_lacosmic(
    gcorr, sigclip=7, readnoise=stdhdr["RDNOISE"] * u.electron, verbose=True
)

Starting 4 L.A.Cosmic iterations
Iteration 1:
325 cosmic pixels this iteration
Iteration 2:
34 cosmic pixels this iteration
Iteration 3:
21 cosmic pixels this iteration
Iteration 4:
19 cosmic pixels this iteration

Show code cell source Hide code cell source

fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
vmin, vmax = interval.get_limits(crrej.data)
im = ax.imshow(crrej.data.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Cosmic Ray Removed Feige34 Image")
fig.show()

Looks good. Now we will identify the source position.

# Slice the image in the spatial direction
lcut = 1000  # left cut (dispersion direction)
rcut = 1050  # right cut (dispersion direction)

apall_probe = np.sum(crrej.data[lcut:rcut, :], axis=0)
x_probe = np.arange(apall_probe.size)
peak_pix = peak_local_max(
    apall_probe, num_peaks=1, min_distance=50, threshold_abs=np.median(apall_probe)
)[0][0]
print("Peak pixel:", peak_pix)

Peak pixel: 68

Show code cell source Hide code cell source

fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(111)
ax.step(x_probe, apall_probe, where="mid")
ax.plot(
    (peak_pix, peak_pix),
    (
        apall_probe[peak_pix] + 0.01 * apall_probe[peak_pix],
        apall_probe[peak_pix] + 0.05 * apall_probe[peak_pix],
    ),
    c="tab:red",
)
ax.annotate(
    peak_pix,
    (peak_pix, apall_probe[peak_pix] + 0.05 * apall_probe[peak_pix]),
    textcoords="offset points",
    xytext=(0, 10),
    ha="left",
    va="top",
    color="k",
    fontsize="small",
    rotation=45,
)
ax.set_title("Aperture Sum")
ax.set_xlabel("Pixel Number (Spatial Direction)")
ax.set_ylabel("Counts [e-]")
fig.show()

Here we crop the image in the dispersion direction, from lcut to rcut. We sum the counts in the dispersion direction to get the 1D spatial profile of the source. We then use the peak_local_max function to identify the source position in the spatial profile. The peak_local_max function finds the local maxima in the 1D profile, which can be used to identify the source position. Then we can also estimate the full width at half maximum (FWHM) of the source, using the peak_widths function from the scipy.signal module.

From the above figure, we were able to get the pixel position of the peak of the source. However, the peak position is not always the center of the source, since it is an integer value. Therefore, we need some centering algorithm to get the center of the source to properly trace the spatial position of the source along the dispersion axis.

Here we use the centroid_com function from the photutils package to get the center of the source. The centroid_com function calculates the center of mass of the source, which is a good approximation of the center of the source. This function originally calculate the centroid of the source in the 2D image, but we will use them to calculate the center of the source in the 1D spatial profile. For more information on centroiding, please refer to the photutils documentation.

Before we estimate the center of the source, we need to estimate the background level. The background level is estimated by fitting a polynomial to the 1D spatial profile of the source. The background level is then subtracted from the 1D spatial profile to get the background-subtracted profile. The background-subtracted profile is then used to estimate the center of the source.

Let’s estimate the background level and subtract it from the 1D spatial profile. Here we use a Chebyshev polynomial of degree 3. The degree of the polynomial can be adjusted depending on the complexity of the background. To model the background, we should set the sky region where the source is not present. The sky region should be chosen carefully to avoid any contamination from the source.

# Select the sky region for the background estimation
sky_limit_peak = 10 # pixel limit from the peak pixel

mask_sky = np.ones_like(apall_probe, dtype=bool)
mask_sky[peak_pix - sky_limit_peak : peak_pix + sky_limit_peak] = False

x_probe_sky = x_probe[mask_sky]
apall_probe_sky = apall_probe[mask_sky]

fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(111)
ax.step(x_probe, apall_probe, where="mid")
ax.fill_between(x_probe, 0, 1, where=mask_sky, alpha=0.1, transform=ax.get_xaxis_transform())
ax.set_xlabel("Pixel Number (Spatial Direction)")
ax.set_ylabel("Counts [e-]")

Text(0, 0.5, 'Counts [e-]')

# Fit the sky background with a Chebyshev polynomial
# define function for background fitting
def fit_background(x, apall, fitmask, sigma_sigclip, order_skymodel):
    x_sky = x[fitmask]
    apall_sky = apall[fitmask]
    # sigma clipping
    mask_sigclip = sigma_clip(apall_sky, sigma=sigma_sigclip, masked=True).mask

    # fitting the sky background
    coeffs_sky, fitfull_sky = chebfit(
        x_sky[~mask_sigclip],
        apall_sky[~mask_sigclip],
        deg=order_skymodel,
        full=True,
    )

    # sky background model
    x = np.arange(apall.size)
    skymodel = chebval(x, coeffs_sky)

    # calculating the rms
    residual_skymodel = fitfull_sky[0][0]
    rms_skymodel = np.sqrt(residual_skymodel / len(x_sky[~mask_sigclip]))

    return skymodel, rms_skymodel, mask_sigclip

# setting the parameters
sigma_sigclip = 3
order_skymodel = 3

# fitting the sky background
skymodel, rms_skymodel, mask_sigclip = fit_background(
    x_probe,
    apall_probe,
    mask_sky,
    sigma_sigclip=sigma_sigclip,
    order_skymodel=order_skymodel,
)

# number of pixels
npix_sky = len(x_probe_sky)
npix_skymodel = np.count_nonzero(~mask_sigclip)

Show code cell source Hide code cell source

fig, axes = plt.subplots(2, 1, figsize=(12, 8), gridspec_kw={"height_ratios": [4, 1]})
ax = axes[0]
ax.step(x_probe, apall_probe, where="mid")
ax.fill_between(
    x_probe, 0, 1, where=mask_sky, alpha=0.1, transform=ax.get_xaxis_transform()
)

ax.plot(
    x_probe,
    skymodel,
    c="tab:red",
    lw=1,
    label=f"Chebyshev Sky Model ({npix_skymodel:d}/{npix_sky:d})",
)
ax.plot(
    x_probe_sky[mask_sigclip],
    apall_probe_sky[mask_sigclip],
    "x",
    c="tab:orange",
    label=f"Sigma-clipped ({sigma_sigclip:.1f}-sigma)",
)

resax = axes[1]
resax.plot([0, x_probe[-1]], [0, 0], c="k", ls="--", lw=0.8)
resax.step(
    x_probe, apall_probe - skymodel, where="mid", label=f"RMS={rms_skymodel:.2f}"
)
resax.fill_between(
    x_probe,
    0,
    1,
    where=mask_sky,
    color="tab:blue",
    alpha=0.1,
    transform=resax.get_xaxis_transform(),
)
resax.plot(
    x_probe_sky[mask_sigclip],
    apall_probe_sky[mask_sigclip] - skymodel[mask_sky][mask_sigclip],
    "x",
    c="tab:orange",
)
resax.set_xlabel("Pixel Number (Spatial Direction)")
resax.set_ylabel("Residuals [e-]")
# ylim_ax = ax.get_ylim()
ylim_resax = rms_skymodel * np.array([-1, 1]) * 10
resax.set_ylim(ylim_resax)
resax.tick_params(axis="x", direction="in", top="on")
resax.legend()

title_str = f"Skyfit: Chebyshev order {order_skymodel:d} ({sigma_sigclip:.1f}-sigma)"
ax.tick_params(axis="x", direction="in", labelbottom=False)
ax.set_xlabel("Pixel Number (Spatial Direction)")
ax.set_ylabel("Counts [e-]")
ax.set_title(title_str)
ax.legend()

for aa in axes:
    aa.set_xlim(0, x_probe[-1])

fig.subplots_adjust(hspace=0)

From the sigma clipping, we can see some rejected points near the source. This indicates that the sky region we set is not sufficiently far from the source. We should set the sky region further away from the source to avoid any contamination from the source. Sigma clipping will also reject any cosmic rays or bad pixels in the sky region, which will affect the background estimation.

from photutils.centroids import centroid_com

apall_probe_skysub = apall_probe - skymodel

center_probe_com = centroid_com(apall_probe_skysub)[0]
print("Center without sky subtraction:", centroid_com(apall_probe)[0])
print("Center after sky subtraction:", center_probe_com)

pwresult = peak_widths(apall_probe_skysub, [peak_pix], rel_height=0.5)
fwhm_probe = pwresult[0][0]
left_ips, right_ips = pwresult[2][0], pwresult[3][0]
fwhm_hight = pwresult[1][0]
print("FWHM:", fwhm_probe)
print("Left and right interpolated positions:", left_ips, right_ips)

Center without sky subtraction: 67.87118517140294
Center after sky subtraction: 67.77667632982664
FWHM: 2.9532888855505917
Left and right interpolated positions: 66.33374812839358 69.28703701394417

Show code cell source Hide code cell source

import plotly.io as pio
import plotly.express as px
import plotly.offline as py
import plotly.graph_objects as go

pio.renderers.default = "notebook"

fig = px.line(
    x=np.arange(len(apall_probe_skysub)),
    y=apall_probe,
    title="Sky Subtracted Sum",
    labels={"x": "Pixel Number (Spatial Direction)", "y": "Counts [e-]"},
    line_shape="hvh"
)

fig.add_vline(
    x=peak_pix,
    line_dash="dash",
    line_color="red",
    annotation_text="Peak",
    annotation_position="top right",
    line_width=0.8,
)

fig.add_vline(
    x=center_probe_com,
    line_dash="dash",
    line_color="green",
    annotation_text="Center",
    annotation_position="top left",
    line_width=0.8,
)

fig.add_trace(go.Scatter(
    x=[left_ips, right_ips],
    y=[fwhm_hight, fwhm_hight],
    mode="lines+markers",
    line=dict(color="red", width=1.5),
    name=f"FWHM={fwhm_probe:.2f}",
))

fig

Zoom in to see the subtle difference between the peak position and the center position.

2. Aperture Tracing

We have identified the source position in the spatial direction by collapsing a slice of the 2D science frame along the dispersion direction. Next, we will trace the spatial position of the source along the dispersion direction. This involves changing the location of the slices in the dispersion direction and identifying the source position in each slice.

# Define the function for determining the center of the source for given slice range
def get_skymask(apall, sky_limit_peak):
    peak_pix = peak_local_max(
        apall,
        num_peaks=1,
        min_distance=10,
        threshold_abs=np.mean(apall),
    )[0][0]

    mask_sky = np.ones_like(apall, dtype=bool)
    mask_sky[peak_pix - sky_limit_peak : peak_pix + sky_limit_peak] = False

    return mask_sky, peak_pix

def find_center(data, lcut, rcut, sky_limit_peak, sigma_sigclip, order_skymodel):
    apall = np.sum(data[lcut:rcut, :], axis=0)
    x = np.arange(apall.size)

    # mask for the sky region and initial peak pixel
    mask_sky, peak_pix = get_skymask(apall, sky_limit_peak)

    # fitting the sky background
    skymodel, rms_skymodel, mask_sigclip = fit_background(
        x,
        apall,
        mask_sky,
        sigma_sigclip=sigma_sigclip,
        order_skymodel=order_skymodel,
    )

    apall_skysub = apall - skymodel

    # center of the source
    center_com = centroid_com(apall_skysub)[0]

    # fwhm of the source

    pwresult = peak_widths(apall_skysub, [round(peak_pix)], rel_height=0.5)
    fwhm = pwresult[0][0]

    return center_com, fwhm


# Find the center of the source for all the slices
slice_width = 50
sky_limit_peak = 20 # increased the distance from the source
midpoints = []
aptrace = []
aptrace_fwhm = []
lcut, rcut = 0, slice_width
while rcut < len(crrej.data):
    rcut = lcut + slice_width
    midpoints.append((lcut + rcut) / 2)
    center, fwhm = find_center(
        crrej.data, lcut, rcut, sky_limit_peak, sigma_sigclip, order_skymodel
    )
    aptrace.append(center)
    aptrace_fwhm.append(fwhm)
    lcut += slice_width
    rcut += slice_width

aptrace = np.array(aptrace)
aptrace_fwhm = np.array(aptrace_fwhm)
aptrace_sigma = aptrace_fwhm * gaussian_fwhm_to_sigma
midpoints = np.array(midpoints)

# Let's fit the aperture trace
order_aptrace = 9
sigma_clip_aptrace = 3

# fitting the trace line
coeffs_aptrace_init = chebfit(midpoints, aptrace, deg=order_aptrace)

# sigma clipping
mask_aptrace = sigma_clip(
    aptrace - chebval(midpoints, coeffs_aptrace_init),
    sigma=sigma_clip_aptrace,
    masked=True,
).mask

# fitting the trace line again
coeffs_aptrace = chebfit(
    midpoints[~mask_aptrace], aptrace[~mask_aptrace], deg=order_aptrace
)
residuals_aptrace = aptrace - chebval(midpoints, coeffs_aptrace)
rejected = ~np.in1d(midpoints, midpoints[~mask_aptrace])
rms_aptrace = np.sqrt(np.mean(residuals_aptrace[~mask_aptrace] ** 2))

# aperture trace model
ximg = np.arange(crrej.data.shape[0])
aptrace_model = chebval(ximg, coeffs_aptrace)

Show code cell source Hide code cell source

# plot the aperture trace
fig, axes = plt.subplots(2, 1, figsize=(12, 6), gridspec_kw={"height_ratios": [3, 1]})
ax = axes[0]
ax.imshow(
    crrej.data.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax, aspect="auto"
)
ax.scatter(midpoints, aptrace, marker="+", c="dodgerblue", s=100, lw=2, label="Data")
ax.scatter(
    midpoints[rejected],
    aptrace[rejected],
    marker="x",
    c="salmon",
    s=100,
    lw=1,
    label="Rejected",
)
ax.plot(ximg, aptrace_model, c="tab:red", lw=1, label="Chebyshev Trace Model")
ax.plot(midpoints, aptrace + 5*aptrace_sigma, c="r", lw=0.8, ls="--", label="5-sigma")
ax.plot(midpoints, aptrace - 5*aptrace_sigma, c="r", lw=0.8, ls="--")
ax.plot(midpoints, aptrace + 10*aptrace_sigma, c="r", lw=0.8, ls=":", label="10-sigma")
ax.plot(midpoints, aptrace - 10*aptrace_sigma, c="r", lw=0.8, ls=":")
ax.set_title("Aperture Trace")
ax.set_ylabel("Pixel Number (Spatial Direction)")
ax.set_xticks([])
ax.legend()

resax = axes[1]
resax.plot([0, ximg[-1]], [0, 0], c="k", ls="--", lw=0.8)
resax.scatter(midpoints, residuals_aptrace, marker="+", c="dodgerblue", s=100, lw=2)
resax.scatter(
    midpoints[rejected],
    residuals_aptrace[rejected],
    marker="x",
    c="salmon",
    s=100,
    lw=1,
)
resax.set_xlabel("Pixel Number (Dispersion Direction)")
resax.set_ylabel("Residuals [pix]")
resax.tick_params(axis="x", direction="in", top="on")
resax_ylim = np.max(np.abs(residuals_aptrace)) * np.array([-1, 1])
resax.set_ylim(resax_ylim)
resax.set_xlim(0, ximg[-1])
resax.text(
    0.99,
    0.95,
    f"RMS={rms_aptrace:.2f}",
    transform=resax.transAxes,
    ha="right",
    va="top",
)

fig.subplots_adjust(hspace=0)

3. Aperture Summation

Now we have traced the spatial position of the source along the dispersion direction. We will sum the flux along the trace to extract the 1D spectrum. The aperture size should be chosen carefully to avoid any contamination from the sky or other sources. The aperture size should be large enough to capture sufficient amount of the flux from the source but small enough to avoid any contamination from the sky or other sources.

See also

The type of aperture used for the extraction can affect the signal-to-noise ratio of the extracted spectrum. One can use a simple aperture, which sums the flux within a fixed aperture size, or an optimal extraction (Horne 1986), which weights the flux based on the spatial profile of the source. The optimal extraction can improve the signal-to-noise ratio of the extracted spectrum, especially when the source is extended or the background is not uniform. In this tutorial, we will use a simple aperture for the extraction, but you can try the optimal extraction for better results.

Let’s make a simple aperture summation to extract the 1D spectrum!

# Define the function for aperture summation
def aperture_sum(
    cut, center, apsum_sigma, ap_sigma, sky_limit_peak, sigma_sigclip, order_skymodel
):
    x = np.arange(len(cut))

    mask_sky, peak_pix = get_skymask(cut, sky_limit_peak)
    # aperture size = trace center +/- apsum_sigma * ap_sigma
    appos_lower = int(center - apsum_sigma * ap_sigma)
    appos_upper = int(center + apsum_sigma * ap_sigma)

    # fitting the sky background
    skymodel, rms_skymodel, mask_sigclip = fit_background(
        x,
        cut,
        mask_sky,
        sigma_sigclip=sigma_sigclip,
        order_skymodel=order_skymodel,
    )

    # subtracting the sky background
    cut_skysub = cut - skymodel
    source = cut_skysub[appos_lower:appos_upper]

    # aperture sum
    ap_sum = np.sum(source)

    # calculating the uncertainty
    sig = rdnoise**2 + source + skymodel[appos_lower:appos_upper]
    ap_sig = np.sqrt(np.sum(sig))

    return ap_sum, ap_sig


# Parameters for the aperture summation
apsum_sigma = 10
ap_fwhm = np.median(aptrace_fwhm[~mask_aptrace])  # median FWHM in pixels
ap_sigma = ap_fwhm * gaussian_fwhm_to_sigma  # sigma under the Gaussian assumption

# Extract the spectrum along the dispersion direction
ap_sum = []
ap_sig = []

for i in range(crrej.data.shape[0]):
    cut_i = crrej.data[i, :]
    center_i = aptrace_model[i]

    ap_sum_i, ap_sig_i = aperture_sum(
        cut_i,
        center_i,
        apsum_sigma,
        ap_sigma,
        sky_limit_peak,
        sigma_sigclip,
        order_skymodel,
    )
    ap_sum.append(ap_sum_i)
    ap_sig.append(ap_sig_i)

ap_sum = np.array(ap_sum) / exptime
ap_sig = np.array(ap_sig) / exptime

Show code cell source Hide code cell source

# Plot the extracted spectrum
fig = go.Figure()
fig.add_trace(go.Scatter(
    x=ximg,
    y=ap_sum,
    error_y=dict(
        type="data",
        array=ap_sig,
        visible=True,
        color="dodgerblue",
        thickness=1,
        width=1,
    ),
    mode="lines+markers",
    marker=dict(color="black", size=3),
    line=dict(color="black", width=0.8),
))
fig.update_layout(
    title="Extracted Spectrum of Feige34",
    xaxis_title="Pixel Number (Dispersion Direction)",
    yaxis_title="Counts per Exposure Time[e-/s]",
    showlegend=False
)
fig

# Let's make a class for the extraction of the spectrum
class InstrumentalSpectrum:
    def __init__(
        self,
        fpath,
        sky_limit_peak,
        sigma_sigclip,
        order_skymodel,
        order_aptrace,
        sigma_clip_aptrace,
        slice_width_aptrace,
        apsum_sigma,
        trim=[0, 4220, 670, 810],
    ):
        self.fname = fpath
        self.TRIM = trim
        self.SKY_LIMIT_PEAK = sky_limit_peak
        self.SIGMA_SIGCLIP = sigma_sigclip
        self.ORDER_SKYMODEL = order_skymodel
        self.ORDER_APTRACE = order_aptrace
        self.SIGMA_CLIP_APTRACE = sigma_clip_aptrace
        self.SLICE_WIDTH_APTRACE = slice_width_aptrace
        self.APSUM_SIGMA = apsum_sigma

        self.hdu = fits.open(fpath)
        self.data = self.hdu[0].data[
            self.TRIM[0] : self.TRIM[1], self.TRIM[2] : self.TRIM[3]
        ]
        self.hdr = self.hdu[0].header
        self.EXPTIME = self.hdr["EXPTIME"]
        self.RDNOISE = self.hdr["RDNOISE"]

        self.ccd = CCDData(data=self.data, header=self.hdr, unit="adu")
        self.gcorr = gain_correct(self.ccd, gain=self.hdr["GAIN"] * u.electron / u.adu)
        self.crrej = cosmicray_lacosmic(
            self.gcorr,
            sigclip=7,
            readnoise=self.hdr["RDNOISE"] * u.electron,
            verbose=False,
        )

        self.set_aptrace()
        self.set_apsum()

    def set_aptrace(self):
        slice_width = self.SLICE_WIDTH_APTRACE
        lcut, rcut = 0, slice_width
        midpoints = []
        aptrace = []
        aptrace_fwhm = []
        while rcut < len(self.crrej.data):
            rcut = lcut + slice_width
            midpoints.append((lcut + rcut) / 2)
            center, fwhm = find_center(
                self.crrej.data,
                lcut,
                rcut,
                self.SKY_LIMIT_PEAK,
                self.SIGMA_SIGCLIP,
                self.ORDER_SKYMODEL,
            )
            aptrace.append(center)
            aptrace_fwhm.append(fwhm)
            lcut += slice_width
            rcut += slice_width

        self.aptrace = np.array(aptrace)
        self.aptrace_fwhm = np.array(aptrace_fwhm)
        midpoints = np.array(midpoints)

        # fitting the trace line
        coeffs_aptrace_init = chebfit(midpoints, self.aptrace, deg=self.ORDER_APTRACE)

        # sigma clipping
        mask_aptrace = sigma_clip(
            self.aptrace - chebval(midpoints, coeffs_aptrace_init),
            sigma=self.SIGMA_CLIP_APTRACE,
            masked=True,
        ).mask
        self.mask_aptrace = mask_aptrace

        # fitting the trace line again
        self.coeffs_aptrace = chebfit(
            midpoints[~mask_aptrace],
            self.aptrace[~mask_aptrace],
            deg=self.ORDER_APTRACE,
        )
        self.residuals_aptrace = self.aptrace - chebval(midpoints, self.coeffs_aptrace)
        self.rejected_aptrace = ~np.in1d(midpoints, midpoints[~mask_aptrace])
        self.rms_aptrace = np.sqrt(np.mean(self.residuals_aptrace[~mask_aptrace] ** 2))

        # aperture trace model
        ximg = np.arange(self.crrej.data.shape[0])
        self.aptrace_model = chebval(ximg, self.coeffs_aptrace)

    def set_apsum(self):
        ap_sum = []
        ap_sig = []
        self.ap_fwhm = np.median(self.aptrace_fwhm[~self.mask_aptrace])
        self.ap_sigma = self.ap_fwhm * gaussian_fwhm_to_sigma

        for i in range(self.crrej.data.shape[0]):
            cut_i = self.crrej.data[i, :]
            center_i = self.aptrace_model[i]

            ap_sum_i, ap_sig_i = aperture_sum(
                cut_i,
                center_i,
                self.APSUM_SIGMA,
                self.ap_sigma,
                self.SKY_LIMIT_PEAK,
                self.SIGMA_SIGCLIP,
                self.ORDER_SKYMODEL,
            )
            ap_sum.append(ap_sum_i)
            ap_sig.append(ap_sig_i)

        self.ap_sum = np.array(ap_sum) / self.EXPTIME
        self.ap_sig = np.array(ap_sig) / self.EXPTIME

std_insspecs = []
for std in std_list_f34:
    fname = OUTDIR / f"p{std.stem}"
    ins = InstrumentalSpectrum(
        fname,
        sky_limit_peak=20,
        sigma_sigclip=3,
        order_skymodel=3,
        order_aptrace=9,
        sigma_clip_aptrace=3,
        slice_width_aptrace=50,
        apsum_sigma=10,
    )
    std_insspecs.append(ins)

Show code cell source Hide code cell source

fig = go.Figure()
f, axes = plt.subplots(3, 1, figsize=(12, 3))
i=0
for StdInsSpec in std_insspecs:
    ax = axes[i]
    fig.add_trace(go.Scatter(
        x=ximg,
        y=StdInsSpec.ap_sum,
        error_y=dict(
            type="data",
            array=StdInsSpec.ap_sig,
            visible=True,
            thickness=1,
            width=1,
        ),
        mode="lines+markers",
        marker=dict(size=3),
        line=dict(width=0.8),
        name=StdInsSpec.fname.stem,
    ))
    i += 1
    ax.imshow(
        StdInsSpec.crrej.data.T,
        cmap="gray",
        origin="lower",
        vmin=vmin,
        vmax=vmax,
    )
    ax.set_title(StdInsSpec.fname.stem)
f.suptitle("2D Images of the Standard Star Frames")
fig.update_layout(
    title="Extracted Spectrum of Feige34",
    xaxis_title="Pixel Number (Dispersion Direction)",
    yaxis_title="Counts per Exposure Time[e-/s]",
    showlegend=True,
    legend=dict(
    yanchor="top",
    y=0.99,
    xanchor="left",
    x=0.01
))
fig

stdtab_f34

Table length=3

filename	date-obs	object	ra	dec	type	exptime	airmass	altitude	disperser	slit-width	slit-pa	slit	filter1	filter2	filter3	ut-str	ut-end
str50	str50	str50	str50	str50	str50	float64	float64	float64	str50	float64	float64	str50	str50	str50	str50	str50	str50
FCSA00141880.fits.gz	2013-06-07	Feige34	10:39:30.972	+43:06:49.26	OBJECT	60.0	1.188	57.2701	SCFCGRHD90	0.5	-240.3	SCFCSLLC08	NONE	SCFCFLSO58	NONE	05:49:08.496	05:50:08.794
FCSA00141894.fits.gz	2013-06-07	Feige34	10:39:30.923	+43:06:49.54	OBJECT	60.0	1.207	55.9261	SCFCGRHD90	0.5	-240.3	SCFCSLLC20	NONE	SCFCFLSO58	NONE	05:58:05.720	05:59:05.978
FCSA00141878.fits.gz	2013-06-07	Feige34	10:39:30.978	+43:06:49.28	OBJECT	60.0	1.185	57.5049	SCFCGRHD90	0.5	-240.3	SCFCSLLC08	NONE	SCFCFLSO58	NONE	05:47:32.380	05:48:32.720

As you can see from the figure above, the extracted spectrum of FCSA00142006 shows different level of continuum. It can be due to various reasons, such as slight shift of the source position with respect to the slit aperture, or the difference in the seeing condition, out-of-focus, etc. To correct for this, we can normalize the extracted spectrum by dividing it by the continuum level at some range. Here we will normalize the extracted spectrum by dividing it by the median value of the continuum level at the range of 1000-2000 pixels.

Note

If we manipulate the continuum level of standard star like this, how should we calibrate the flux of the science frame and trust them? For this reason, we need a process called aperture correction to correct the flux. We will skip this process for simplicity, but you should always consider this process in your own data reduction.